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METHOD AND SYSTEM FOR SOLVING FINITE ELEMENT MODELS 
USING MULTI-PHASE PHYSICS 

Related Applications 

[0001] This application claims priority under 35 U.S.C. 
§ 119(e) to United States Patent Application 
No. 60/215,697 entitled "Method and System for Oil 
Reservoir Simulation and Modeling" by Stephen R. Kennon, 
Kok Thye Lim, Scott A. Canaan, Steven B. Ward, Stuart W. 
Pond, Jr. and Edward J. Barragy, filed June 29, 2 0 00, 
which is incorporated by reference as if set forth in its 
entirety herein. 

Field of the Invention • 

[0002] The invention relates generally to methods for 

modeling physical systems using finite element analysis 
and, more specifically, to methods for solving finite 
element models wherein the matrix that controls the 
solution for the model is configured to maintain the 
properties of monotonicity and linearity preservation. 

Background of the Invention. 

[0003] Physical systems can be modeled mathematically to 

simulate their behavior under certain conditions. There 
are a wide variety of means to model these systems, 
ranging from the very simplistic to the extremely 
complicated. One of the more sophisticated means to model 
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physical systems is through the use of finite element 
analysis. As the name implies, finite element analysis 
involves the representation of individual, finite elements 
of the physical system in a mathematical model and the 
solution of this model in the presence of a predetermined 
set of boundary conditions. Another comparable means to 
model physical systems is through the use of finite 
difference analysis. Finite difference analysis involves 
the modeling of individual points within a modeled space 
and computing the differences between these points., 
Finite difference analysis is often used for simulating 
the dynamic behavior of fluids. 

[0004] Traditional finite difference techniques and finite 
element techniques using streamlined, upwinding methods 
are typically used to simulate the production of oil in a 
reservoir. While each of these techniques has its own 
advantages, it also has its own disadvantages. Generally 
speaking, the finite difference techniques produce 
physically realistic results, but they are not very 
accurate. The finite element techniques, on the other 
hand, are more accurate, but they produce results which 
are not physically realistic. 

[0005] As a result of their respective disadvantages, both 
the finite difference techniques and finite element 
techniques which are conventionally used normally require 
a great deal of computer resources. In the case of finite 
difference techniques, reasonable accuracy can be 
achieved, but this requires many more nodes then would be 
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necessary in a finite element model. This increases the 
amount of memory and CPU time which are needed to compute 
an accurate solution. In the case of finite element 
techniques, comparable accuracy can be achieved with fewer 
nodes, but the solutions may not be realistic. For 
example, concentrations may be more than 100 percent, or 
permeabilities may be negative. Consequently, there may 
be problems for which the solution does not converge, or 
for which the solution may converge very slowly. These 
techniques are therefore less reliable and may require a 
large number of iterations before acceptable accuracy is 
achieved. 

[0006] Because each of these standard techniques has its own 
drawbacks, and because these drawbacks increase in the 
amount of computer resources which are necessary to 
generate acceptable solutions, it would be desirable to 
provide a method for modeling systems such as oil 
reservoirs which reliably produces accurate, realistic 
solutions for these systems. 

Summary of the Invention 

[0007] One or more of the problems outlined above may be 
solved by the various embodiments of the invention. 
Broadly speaking, the invention comprises systems and 
methods for solving finite element models, wherein the 
matrix that governs the solution is modified so that it 
satisfies the properties of monotonicity and preservation 
of linearity. This is accomplished by adjusting the 
weighting coefficients of the matrix so that the elements 
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which lie on the diagonal of the matrix are non-negative 
and the elements which are off the diagonal are non- 
positive. This causes the solution to be both more 
accurate than traditional finite difference techniques and 
more realistic than traditional finite element techniques, 

[0008] The present techniques may provide various advantages 
over the prior art. For instance, because the present 
techniques are more accurate, the number of iterations 
which are required to produce an acceptable result may be 
reduced, thereby conserving computing resources. The 
present techniques may also cause solutions to more 
reliably converge, since they are constrained to the range 
of realistic results and do not oscillate in the same 
manner as the solutions of traditional finite element 
techniques. The present techniques are therefore 
generally faster, more accurate, less resource -intensive 
and more reliable than prior art techniques. 

[0009] One embodiment of the present invention comprises a 
method for generating a solution for a modeled system 
using finite element techniques, wherein the system is 
discretized on a finite element mesh and wherein the 
contribution of each node to the discretization is 
weighted based upon the direction of fluid flow across 
each element. The contributions are weighted to favor the 
nodes which are upstream from the other nodes of the 
respective elements. Weighting the contributions of the 
nodes in this manner causes the resulting matrix to be 
essentially diagonal. That is, the elements which lie on 
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the diagonal of the matrix are emphasized, while the 
elements which are off the diagonal are de- emphasized. 
Once the matrix has been modified in this manner, it is 
solved using traditional techniques. 

[0010] An alternate embodiment of the present invention 

comprises a method for modeling the flow of multiphase 
fluids in an oil reservoir. In this embodiment, a four 
dimensional finite element mesh corresponding to the oil 
reservoir is constructed. A solution is generated for the 
mesh by distributing the residuals of a finite element 
operator among the nodes of the mesh and then assembling 
and solving the system using traditional finite element 
techniques. In this embodiment, distributing the residual 
of the finite element operator comprises, for each 
element, computing the value of the finite element 
operator for each node, dividing these values into a set 
of positive values and a set of negative values, 
distributing the residual to the nodes to emphasize the 
diagonal elements and de-emphasize the off-diagonal 
elements and then computing a tangent operator for the 
matrix. 

[0011] Another alternate embodiment of the present invention 
comprises a software product which is embodied in a 
computer readable medium. The software product comprises 
a plurality of instructions which are configured to cause 
a computer or any other data processor to perform the 
steps of a method in accordance with the present 
disclosure. The computer readable medium may comprise any 
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of a number of different media, including but not limited 
to CD ROMs, DVDs, floppy disks, computer memories, 
magnetic tapes, etc. Still another alternate embodiment 
of the present invention comprises a computer system which 
is configured to perform the steps of a method in 
accordance with the present disclosure. The computer 
system may be programmed to perform such a method by 
instructions which are contained in hardware, firmware or 
software . 

[0012] Numerous additional alternative embodiments are also 
possible . 

Brief Description of the Drawings 

[0013] Other objects and advantages of the invention may 
become apparent upon reading the following detailed 
description and upon reference to the accompanying 
drawings . 

[0014] FIGURE 1 is a diagram illustrating simplices in two 
and three dimensions . 

[0015] FIGURE 2 is a diagram illustrating the difference 

between the actual values in the physical system and the 
solutions generated using traditional finite element 
techniques and techniques according to the present 
invention. 

[0016] FIGURE 3 is a diagram illustrating the flow of fluid 
across a first two-dimensional simplex element. 
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[0017] FIGURE 4 is a diagram illustrating the flow of fluid 
across the first two-dimensional simplex element in a 
direction opposite the flow illustrated in FIGURE 3 . 

[0018] FIGURE 5 is a diagram illustrating the flow of fluid 
across a second two-dimensional simplex element. 

[0019] FIGURE 6 is a diagram illustrating the flow of fluid 
across the second two-dimensional simplex element in a 
direction opposite the flow illustrated in FIGURE 5. 

[0020] FIGURE 7 is a diagram illustrating a two-dimensional 
projection of a time-space domain. 

[0021] FIGURE 8 is a flow diagram illustrating a method in 
accordance with one embodiment of the present invention. 

[0022] While the invention is subject to various 

modifications and alternative forms, specific embodiments 
thereof are shown by way of example in the drawings and 
the accompanying detailed description. It should be 
understood, however, that the drawings and detailed 
description are not intended to limit the invention to the 
particular embodiment which is described. This disclosure 
is instead intended to cover all modifications, 
equivalents and alternatives falling within the scope of 
the present invention as defined by the appended claims. 

Detailed Description of a Preferred Embodiment 

[0023] A preferred embodiment of the invention is described 
below. It should be noted that this and any other 
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embodiments described below are exemplary and are intended 
to be illustrative of the invention rather than limiting. 

[0024] Broadly speaking, the invention comprises methods for 
solving finite element models, wherein the solution is 
guaranteed both to be monotonic and to satisfy linearity 
preservation. This ensures that the solution will be both 
accurate and realistic. The monotonicity and the 
linearity of the solution are ensured by generating a 
matrix for the solution of the finite element model such 
that the elements of the matrix which are on the diagonal 
are non-negative, and the elements which are off -diagonal 
are non-positive. This matrix can then be solved using 
conventional finite element techniques to generate the 
desired accurate, realistic solution. 

[0025] The present methods provide a number of advantages 
over traditional finite element and finite difference 
techniques. As indicated above, the present methods 
provide the accuracy of traditional finite element 
techniques without the disadvantages of unrealistic, 
oscillating solutions. Similarly, the present methods 
provide the stability of traditional finite difference 
techniques while providing much greater accuracy. 

[0026] The accuracy and realism provided by the present 

methods lead to additional advantages in the utilization 
of computer resources. For example, the present methods 
may allow for faster generation of solutions for physical 
systems. This is a result of several things. First, 
because the solution is physically realistic, it does not 
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oscillate, and consequently tends to converge more 
rapidly. Further, because generation of the solution 
using the present methods typically requires less 
iterations and is faster than traditional finite element 
techniques, less CPU time is required to generate the 
solution. The present methods may also use less memory 
than traditional techniques. Because the present methods 
are generally more accurate than traditional finite 
difference techniques, is not necessary to use as many 
nodes to model the physical system. Consequently, less 
memory is required to store the model (and less CPU time 
is required to compute a solution based on these nodes.) 
the present methods may also be more reliable than 
traditional techniques. The oscillation of solutions 
generated using traditional finite element techniques may, 
in some instances, prevent a solution from converging at 
all. As noted above, the present methods provide 
solutions which are realistic (e.g., they do not oscillate 
between realistic solutions and solutions which are not 
possible,) so the solutions will reliably converge. 

[0027] In finite element modeling, the region that is to be 
analyzed is broken up into sub-regions called elements. 
This process of dividing the region into sub-regions may 
be referred to as discretization or mesh generation. The 
region is represented by functions defined over each 
element. This generates a number of local functions that 
are much simpler than those which would be required to 
represent the entire region. The next step is to analyze 
the response for each element. This is accomplished by 



Gray Cary\AU\4 062731 . 1 
104335-990000 



ATTORNEY DOCKET NO. : PATENT APPLICATION 

OBJlllO-8 CUSTOMER ID: 25094 

10 



building a matrix that defines the properties of the 
various elements within the region and a vector that 
defines the residuals acting on each element in the 
domain. Once all the element matrices and vectors have 
been created, they are combined into a matrix equation. 
After applying boundary conditions, the matrix equation 
can be solved to obtain unknown nodal responses. Intra- 
element responses can be interpolated from nodal values 
using the functions which were defined over each element. 

[002 8] The following description of one embodiment of the 

present method relates to the modeling of oil reservoirs. 
This embodiment is intended to be exemplary, and it should 
be noted that the present methods may be applied to a 
variety of physical systems, such as flow over aircraft, 
weather prediction, flow through internal combustion 
engines, and other industrial fluid flows. This list of 
potential applications of the present methods is intended 
to be illustrative rather than limiting. 

[0029] The details of a preferred embodiment will be set 

forth below. It may be helpful, however, to first define 
a few terms. 

[0030] A node is a point in space. In finite element 

modeling, nodes form the vertices of the elements which 
are modeled. The nodes also form part of a mesh of nodes 
and edges which define the boundaries between elements in 
the modeled space. 
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[0031] An edge is a line between two nodes which form 

vertices of an element. The edges form part of the mesh 
which defines the boundaries between elements in the 
modeled space. 

[0032] A simplex is a spatial configuration of n dimensions 
determined by n + 1 points in a space of dimension equal 
to or greater than n. In other words, a-simplex is a 
geometric spatial element having the minimum number of 
boundary points necessary to enclose a space in a given 
number of dimensions. For example, in two dimensions, a 
simplex comprises a triangle, together with the interior 
area bounded by the triangle (see FIGURE 1.) Two points 
are insufficient to form a simplex in two-dimensional 
space because no area is bounded by the points (and the 
lines which interconnect them.) While four points may be 
sufficient to bound a two-dimensional area, they do not 
comprise the minimum number of boundary points by which 
the two-dimensional area can be bounded. In three 
dimensions, a simplex comprises a tetrahedron, which is 
bounded by four vertices (see FIGURE 1.) In four 
dimensions, a simplex comprises a hypertetrahedron 
(sometimes referred to as a hypertet) having five 
vertices . 

[0033] A mesh is a collection of elements that fill a space. 
These elements are representative of a system which 
resides in that space. Because each element can be 
defined by a plurality of nodes and/or the edges between 
those nodes, a mesh may alternatively be considered a 
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collection of nodes and/or the edges between them. At 
various points in this disclosure, "mesh" will be used to 
alternately refer to collections of elements or 
nodes/edges, depending upon the context in which the term 
is used. The mesh may also be referred to herein as a 
finite element model or simply a model. 

[0034] There is a fundamental difference between modeling 

systems in which there is a single fluid phase and those 
in which there are multiple phases. When there are 
multiple fluid phases, there is a convection and 
potentially a mixing of the phases. (It is assumed that 
there is a pressure which drives the fluids and causes the 
front between them to be transported.) The finite element 
model must therefore include convective operators or 
components. These operators would not be necessary if 
only one fluid phase were modeled, or if the multiple 
fluid phases were completely mixed or emulsified. 

[0035] In a system having multiple fluid phases, 

characteristics of the system may change very abruptly. 
This situation is illustrated in FIGURE 2. FIGURE 2 is a 
diagram illustrating the difference between the actual 
values in the physical system and the solutions generated 
using traditional finite element techniques and techniques 
according to the present invention. It can be seen that 
the actual ("real world") values 11 form what is 
essentially a step function, changing from a first value 
vl to a second value vO immediately upon reaching a 
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distance dl. In this example, the values vO and vl define 
the limits of the range of physically realizable values. 

[0036] When traditional finite element techniques are used to 
model such a system, the solution will typically tend to 
oscillate. It can be seen in FIGURE 2 that the solution 
12 generated using these techniques oscillates above and 
below the actual physical values which it is attempting to 
model. Two of the most significant problems with this 
solution are the obvious inaccuracy (the solution does not 
closely track the actual values) and the fact that, at 
some points, the solution is greater than the maximum 
physically realizable value or less than the minimum 
physically realizable value. The oscillations which occur 
in solutions generated using traditional finite element 
techniques may therefore produce unrealistic results, or 
may even fail to converge at all. 

[0037] Methods in accordance with this disclosure, however, 
do not suffer from these disadvantages because the 
solutions are constrained to be realistic. This 
constraint is achieved by ensuring that the matrix which 
governs the solution vector (i.e., the Jacobian matrix) 
satisfies the properties of monotonicity and preservation 
of linearity. In more practical terms, this consists of 
ensuring that all of the coefficients on the diagonal of 
the matrix are non-negative and that all of the 
coefficients off the diagonal of the matrix are non- 
positive. When the matrix has these characteristics, the 
solution remains realistic at all points, and more 



Gray Cary\AU\4 062731 . 1 
104335-990000 



ATTORNEY DOCKET NO. 
OBJlllO-8 



14 



PATENT APPLICATION 
CUSTOMER ID: 25094 



accurately tracks the behavior of the physical system, as 
indicated by curve 13 in FIGURE 2. The model is therefore 
more likely to produce accurate results and to reliably 
converge . 

[003 8] The Jacobian matrix is derived from the partial 
differential equations that describe multiphase flow 
through porous media based on Darcy' s law, and 
representing conservation of mass and momentum in the 
system. These equations are discretized onto the finite 
element mesh, and the coefficients of the equations form a 
matrix. Each coefficient will be either positive or 
negative, depending upon whether the gradient of the 
equation is positive or negative (i.e., whether the fluid 
flow is positive or negative) . If the finite element 
technique is applied to the differential equations, a 
large set of algebraic equations is produced. That set of 
algebraic equations is nonlinear. The set of equations 
can be solved by using Newton's method (which involves 
taking the first derivative of the equations) . The first 
derivative is the Jacobian matrix. The coefficients of 
this matrix are weighted to obtain a matrix which has non- 
negative terms on its diagonal and non-positive terms off 
the diagonal. The appropriate weighting of the terms 
ensures the monotonic and linear properties of the matrix. 

[0039] The coefficients of the Jacobian matrix are 

essentially weighting factors for each of the nodes of the 
elements in the finite element mesh. In traditional 
finite element techniques, the nodes of an element are 
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evenly weighted. For example, in a two-dimensional 
element, each of the three nodes would be weighted by 1/3. 
If a three-dimensional element were used, each of the four 
nodes of the element would be weighted by 1/4. In the 
present methods, however, the nodes are weighted to 
reflect the flow of the fluids in the modeled system. 
This is illustrated with reference to FIGURES 3-6. 

[0040] FIGURE 3 shows a single, two-dimensional element (a 

triangle) . The flow of fluid with respect to the element 
is indicated by the arrows, which show that the fluid is 
flowing from the upper, right-hand portion of the figure 
to the lower, left-hand portion of the figure. Using 
traditional finite element techniques, each of nodes 21 - 
23 would be evenly weighted (i.e., 1/3, 1/3, 1/3). Using 
the present techniques, the nodes are not evenly weighted. 
Because the fluid is flowing generally from nodes 22 and 
23 to node 21, each of nodes 22 and 23 is weighted by 1/2, 
while node 21 is weighted by 0. Referring to FIGURE 4, it 
can be seen that the weighting of the nodes is different 
when the direction of the fluid flow changes. In FIGURE 
4, the fluid flow is generally from the lower, left-hand 
portion of figure to the upper, right-hand portion of the 
figure. Consequently, node 21 is weighted by 1, while 
nodes 22 and 23 are weighted by 0. 

[0041] Referring to FIGURES 5 and 6, it can be seen that the 
weighting of the nodes is also affected by the shape of 
the corresponding element. In FIGURE 5, the flow of fluid 
is approximately perpendicular to an edge of the element. 
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This situation is similar to that of FIGURE 3, in which 
each of the nodes which terminates the forward edge of the 
element was weighted by 1/2. In that case, the remaining 
node of the element was downstream from a midway point on 
the forward edge of the element. In FIGURE 5, the 
remaining node is instead directly downstream from one of 
the first to nodes. In this situation, the node 32 which 
is upstream from the last node 33 is weighted by 2/3, 
while the other upstream node 31is weighted by 1/3. In 
FIGURE 6, there is a single upstream node 33. Nodes 31 
and 32 are both on the trailing edge of the element, so 
they are both weighted by 0. The upstream node 33, the 
other hand, is weighted by 1. 

[0042] In each of these cases, the weighting of the nodes 
with respect to the direction of fluid flow effectively 
emphasizes the diagonal elements of the matrix, while 
reducing or eliminating the effect of the off -diagonal 
elements. This, in turn, ensures that the matrix has the 
properties of monotonicity and linearity. 

[0043] Although the foregoing is sufficient to enable a 

person of ordinary skill in the field of finite element 
analysis to implement the present techniques, a brief 
technical discussion of the techniques follows. 

[0044] Considering a physical system which is to be modeled 
over time, a finite element mesh is first constructed 
(i.e., the system is discretized in the space-time 
domain) . Referring to FIGURE 7, a two-dimensional 
projection of the time-space domain is illustrated. 
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[0045] ^ ^pace^t^'T] space-time domain to be modeled, 

Q £91^ 

where 

[0046] ^ is broken down (segmented) into a plurality of non- 
overlapping sub-domains , Eacli of these sub-domains 

and e comprises a four- dimensional 

simplex (a hypertetrahedron) . A linear finite element 
basis is formed on each element (^e) , and trial functions 
and test functions are defined on each element. 

[0047] The flow of multi-phase fluids in porous media is 
governed by the equation: 
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Where: 



i = 


1,2, ...,nphases 




= -K(VPi-hPigVZ) 


p 


= P 




= P-P^,, i=2, ...,nphases 












= Viscosity of i'th phase 


Bj 


" Formation Volume Factor of i'th phase 


Pi 


- Density of i'th phase 




= Rock Porosity 


Z 


= Negative of Depth 


Si 


= Saturation of i'th phase 



nphase 



S Sj = 1 

i 

[0048] 



[0049] A standard GalerJcin Finite Element discretization of 



this multi-phase flow equation gives, for each element; 



[0050] 
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[0051] The basic method of employing the present techniques 
is defined by the way that '^i is defined. The present 
method is therefore defined as follows: 

[0052] 

[0053] Referring to FIGURE 8, a flow diagram illustrating a 
method in accordance with one embodiment of the present 
invention is shown. As indicated above, when the 
differential equations are discretized on the finite 
elements, they are basically converted from differential 
equations to algebraic relationships between the values of 
the nodes. This is accomplished by assuming that there is 
a linear variation of the corresponding property (e.g. 
saturation, pressure, permeability, prostate, etc.) 
between the nodes. This restricts the space of all 
possible solutions to only those solutions which are 
linear. When this assumption is used, a set of relations 
are derived from the differential equations for values of 
the coefficients which are linear over each element. When 
this is done, each element contributes a certain amount of 
the discretization to each one of its nodes. The 
contributions which are summed at each node are the 
residuals. It is desired to make the residuals equal 0 in 
order to solve the nonlinear algebraic system. 
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[0054] To distribute the residuals of the finite element 

operator ( J) / Kj is first computed for each node of the 4- 
tet (hypertetrahedron) . The resulting Kj are divided into 

Y- 

positive (Kj^) and negative (Kj") values. The residual J 
is then distributed to the nodes, and the tangent operator 
(Tjk) is computed. The system is then assembled and solved 
in the same manner as if an ordinary finite element 
technique had been used. (Since the process by which 
finite element systems are assembled and solved is well 
known, it will not be described in further detail here) . 

[0055] The benefits and advantages which may be provided by 
the present invention have been described above with 
regard to specific embodiments. These benefits and 
advantages, and any elements or limitations that may cause 
them to occur or to become more pronounced are not to be 
construed as a critical, required, or essential features 
of any or all of the claims. As used herein, the terms 
"comprises," "comprising," or any other variations 
thereof, are intended to be interpreted as non-exclusively 
including the elements or limitations which follow those 
terms. Accordingly, a process, method, article, or 
apparatus that comprises a list of elements does not 
include only those elements but may include other elements 
not expressly listed or inherent to the claimed process, 
method, article, or apparatus. 

[0056] While the present invention has been described with 
reference to particular embodiments, it should be 
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understood that the embodiments are illustrative and that 
the scope of the invention is not limited to these 
embodiments. Many variations, modifications, additions 
and improvements to the embodiments described above are 
possible. Particularly, these variations may include 
computers or other data processing devices, computer 
readable media (such as floppy disks, CD-ROMs, DVD-ROMs, 
etc.,) storage devices, computer memories and the like 
which contain software, firmware or other programming 
embodying the foregoing methods. It is contemplated that 
these variations, modifications, additions and 
improvements fall within the scope of the invention as 
detailed within the following claims. 
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